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MODIFICATIONS  TO  SOLA-VOF  FOR  FLOW  DYNAMICS  IN  SPINNING  CYLINDERS 


1.  INTRODUCTION 

/ 

^->F1uid  dynamic  processes  within  spinning  containers  involve  many 
complex  phenomena.  Boundary  layer  effects  and  free  surface  motions  can  be 
strongly  coupled  with  angular  momentum  effects  to  produce  highly  complex  flow 
structures.  This  situation  exists  in  fluid-filled  projectiles.  -Experimental 
studies  indicate  that  a  better  predictive  capability  is  needed  to  resolve  a 
host  of  important  questions  -  for  instance,  how  fluid  motions  affect  projectile 
stability,  how  spinning  fluid  is  ejected  from  a  projectile  with  an  open  end, 
and  how  two  fluids  initially  separated  within  a  projectile  will  mix  when  the 
separation  device  is  removed. 


The  beginnings  of  a  powerful  analysis  tool  that  can  answer  these 
questions  is  available  in  the  SOLA-VOF  computer  program. This  program  numeri¬ 
cally  solves  the  Navier-Stokes  equations  for  an  incompressible  fluid.  Its 
principal  advantage  over  other  programs  is  its  capability  for  describing  the 
dynamics  of  fluids  with  complicated  free  boundary  behavior.  Unfortunately, 
this  code  does  not  have  provisions  for  an  azimuthal  (swirl)  velocity  component. 
To  remedy  this,  the  project  that  is  described  in  this  report  was  undertaken  to 
extend  SOLA-VOF  to  have  an  azimuthal  velocity  component  that  could  be  used  iri 
either  the  one-fluid  with  free  surface  or  two-fluid  with  interface  cases.  The 
azimuthal  velocity  addition  is  also  to  be  compatible  with  all  other  capabilities 
of  the  original  code. 


For  efficiency,  in  some  calculations  we  found  it  necessary  to  use  a 
new  algorithm  called  SADI  for  pressure  solution.  This  algorithm  is  a  variant  of 
the  well-known  Alternating  Direction  Implicit  (ADI)  method.  SADI  is  more 
versatile  than  ADI,  however,  because  its  implicitness  can  be  limited  to  just 
the  x-direction  or  just  the  y-direction.  In  any  case,  the  method  generally 
converges  significantly  more  rapidly  than  the  standard  ADI  method. 


In  the  course  of  this  project,  we  have  also  included  in  the  new  code, 
SOLA-VOF/CSL ,  a  variety  cf  miscellaneous  updates  and  corrections  to  the  original 
SOLA-VOF  program. 


In  the  next  section  the  theory  for  incorporating  an  azimuthal  velocity 
component  into  SOLA-VOF  is  briefly  described.  Then,  in  Section  3,  a  discussion 
is  given  of  the  numerical  implementation  of  this  theory.  Section  4  contains  a 
description  of  the  new  input  parameters  for  the  code  that  ai  e  required  for 
problems  involving  rotating  fluids.  Finally,  in  Section  5,  two  sample  problems 
are  presented  that  illustrate  the  added  capabilities. 


^Nichols,  B.  D. ,  Hirt,  C.  W. ,  and  Hotchkiss,  R.  S.  SOLA-VOF:  A  Solution 
Algorithm  for  Transient  Fluid  Flow  with  Multiple  Free  Boundaries.  Los 
Alamos  Scientific  Laboratory  Report  LA-8355.  1980. 

^Hirt,  C.  W. ,  and  Nichols,  B.  D.  Volume  of  Fluid  (VCF)  Method  for  the  Dynamics 
of  Free  Boundaries.  J.  Computational  Physics  ^9,  201  (1981). 


2. 


THEORY 


The  original  SOLA-VOF  code  was  primarily  designed  for  incompressible 
flow  problems  governed  by  the  Navler-Stokes  equations 
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Two  dimensional  flow  is  assumed  in  either  plane  geometry  (x,y)  or  cylindrical 
geometry  (r,z).  The  parameter  value  5  -  0  is  the  setting  for  the  plane  case, 
while  c  =  1  is  the  setting  for  cylindrical  coordinates  (i.e. ,  axisymmetric 
flow). 


A  feature  of  SOLA-VOF  that  sets  it  apart  from  other  incompressible  flow 
codes  is  its  ability  to  treat  complex  free  surface  motions  (or  two-fluid  inter¬ 
faces).  This  code  also  contains  many  user  convenient  features  and  has  been 
extensively  utilized  for  a  wide  range  of  problems.  These  aspects  of  SOLA-VOF 
make  it  a  highly  useful  and  reliable  tool  for  the  analysis  of  complicated  flow 
problems. 


To  extend  the  code  for  rotating  (or  swirling)  flow,  we  must  add  an 
additional  momentum  equation  governing  the  azimuthal  velocity,  w,  and  a  centrifu¬ 
gal  acceleration  term  must  be  added  to  the  u-equation,  Eq.  (1).  Still  assuming 
axisymmetric  flow  conditions,  the  new  u-  and  w-equations  are:^ 
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^Landau,  L.  D. ,  and  Lifshitz,  E.  M.  Fluid  Dynamics.  Pergamon  Press,  London. 
1959. 


Here  we  have  retained  the  5  multiplier  on  the  uniquely  cylindrical  terms,  even 
though  the  w-equation  has  no  meaning  for  azimuthal  velocity  in  the  plane  geometry 
case,  5=0.  The  reason  for  this  is  to  maintain  flexibility  in  the  basic  code 
so  that,  for  example,  the  w-equation  could  still  be  used  for  temperature 
or  some  other  scalar  quantity  when  the  code  is  used  for  plane  coordinate  calcula¬ 
tions. 

The  axial  momentum  equation  and  the  equation  for  incompressibility 
are  unchanged  in  the  presence  of  an  azimuthal  velocity. 

Boundary  conditions  on  the  azimuthal  velocity  are  w  =  wg  at  a  no-slip 
wall,  where  wg  is  the  azimuthal  wall  velocity.  At  a  free-slip  wall  there  is  a 
zero  shear  stress,  which  implies  that  3w/3z  =  0  along  a  wall  normal  to  the  z-axis 
while  along  a  cylindrical  wall  (r  =  constant)  we  must  have  3(w/r)/ar  =  0.  The 
code  also  allows  for  continuative  boundaries,  constant  pressure  boundaries,  and 
periodic  boundaries.  For  continuative  and  constant  pressure  boundaries,  the 
free-slip  wall  conditions  are  used  for  the  azimuthal  velocity  as  this  seems  to 
be  the  most  reasonable  choice  for  these  cases. 

When  free  surfaces  are  present,  we  have  simply  assumed  that  3w/an  =  0, 
where  the  derivative  is  in  the  direction  normal  to  the  surface.  For  very  low 
Reynolds  number  flows,  say  Re  <  10,  all  the  free  surface  stress  conditions  used 
in  S0LA-V0F  should  be  revised. ^ 

3.  NUMERICAL  CONSIDERATIONS 

3 . 1  Azimuthal  Momentum  Equation. 

Numerical  approximations  of  Eq.  (4)  must  be  made  with  some  care.  As 
noted  in  the  original  description  of  S0LA-V0F,1  we  cannot  use  a  strictly  conserva 
tive  formulation  of  the  advection  terms  because  in  a  nonuniform  mesh  our  differ¬ 
ence  approximations  would  be  zeroth  order  accurate.  However,  we  want  to  have 
as  close  an  approximation  to  the  conservative  case  as  possible.  The  following 
discussion  will  illustrate  the  approach  we  have  taken. 

Consicer  the  simple  advection  problem  in  conservation  form 

3Q/3t  +  3uQ/3x  +  3v(j/3y  =  0,  (5) 

where  Q  is  any  scalar  quantity.  We  assume  a  staggered  mesh  in  which  Q  values  are 
located  at  center?  of  mesh  cells  and  velocities  are  located  at  cell  edges 
(Figure  1). 


^Hirt,  C.  W. ,  and  Shannon,  J.  P.  Free  Surface  Stress  Conditions  for  Incompress¬ 
ible  Flow  Calculations.  J.  Computational  5hysics  2,  403  (1968). 

^Nichols,  B.  D. ,  and  Hirt,  C.  W. ,  Improved  Free  SurTace  Boundary  Conditions 
for  Numerical  Incompressible  Flew  Calculations.  J.  Compuutational  Physics 
8,  434  (1971). 
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Q. 

Figure  1.  Staggered  Mesh 

Equation  (5)  is  written  in  conservative  form.  The  corresponding  non- 
conservative  advection  equation  is 

3Q/3t  +  u(3Q/3x)  +  v(3Q/3y)  =  0,  (6) 

where  we  have  assumed  the  fluid  Is  incompressible, 

3u/ 3x  +  3v/ 3y  =  0. 

Our  objective  is  to  write  a  finite-difference  approximation  to  Eq. (6) 
that  retains  at  least  first-order  accuracy  In  a  nonuniform  mesh.  Furthermore, 
we  would  like  this  approximation  to  reduce  to  one  for  Eq.  (5)  when  the  mesh  has 
uniform  cell  sizes.  To  do  this  we  begin  Dy  writing  an  approximation  to  the 
conservative  equation  (Eq.  G) 

(QCn  "  QS)/6t  +  (uRQDR  “  uLQDL^*xC  +  (vTQDT  “  vBQDB^6yC  =  0  * 

where  superscripts  refer  to  time  levels  and  6xc  and  6yc  are  the  width  and  height 
of  the  cell  (see  Fig.  1  for  definitions  of  the  subscript  notation).  Here  the  Q 
values  chosen  for  fluxing  across  each  cell  boundary  are  weighted  between  centered 
and  donor  (or  upstream)  values,  e.g. , 

QDR  =  CQc  +  qR  ■  asgn(uRHQR  -  Qc)]/2  (8) 


The  case  a  *  0  gives  a  centered  approximation  and  a  =  1  the  donor  cell  approxima 
tion.  The  expression  sgn(u^)  means  the  sign  of  ur. 

The  continuity  equation  for  each  cell  is 

(UR  -  ul)/5xc  +  (vT  -  vB)/$yc  =  0  (9) 

By  subtracting  Qc  times  Eq.  (9)  from  Eq.  (7)  we  find 

(Q£  1  -  Qj)/«5t  +  ur(Qdr  -  Qc)/5xc  +  ul(Qc  -  Qdl)/6xc  (10) 

+  vt(Qdt-  Qc)/5yc  +  vb(qc-  QDB)/6yc  =  0 
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This  expression  is  an  approximation  to  Eq.  (6).  Unfortunately,  when  the  mesh 
is  nonunifonn  and  we  select  an  a  value  different  from  zero  (which  we  must  for 
numerical  stability  reasons)  the  approximations  in  Eq.  (7),  or  equivalently  £q. 
(10),  are  zeroth  order  accurate.  The  reason  for  this  is  that  the  various  6x 
and  6y  contributions  in  a  Taylor  series  expansion  will  combine  in  such  a  way  as 
to  leave  coefficients  on  the  zeroth  order  advective  terms  that  are  different 
from  unity.  For  instance,  the  term  u3Q/3x  will,  for  pure  donor  cell 
differencing  (a  =  1  and  u  >  0),  become 


u(3Q/3x)(6xl  +  6xc)/6xc 

The  ratio  of  cell  sizes  is  generally  zeroth  order  (it  can  be  first  order  in 
special  circumstances,  depending  on  the  construction  of  the  mesh). 

To  eliminate  this  problem  we  have  replaced  Eq.  (10)  by 


(Qc  Qc)/5t  +  ur(Qdr  Qc^5xRC  +  Ul/3C  "  QDl//(5xLC 


+  vt(qdt  -  Qc)/5yTC  +  vb(Qc  -  QDB)/5yBC  =  0 


(11) 


where,  e.g. ,  fixp^  =  (jxp  +  sxq)/2.  With  this  small  change  in  the  «x  and  sy 
divisions  we  no  longer  have  a  strictly  conservative  approximation  (unless  the 
mesh  is  uniform),  but  we  retain  first  order  spatial  accuracy  even  when  a  *  0. 

Thus,  Eq.  (11)  is  the  recommended  finite-difference  approximation  for  the  azimuthal 
velocity  equation,  which  has  been  used  in  S0LA-V0F/CSL. 


Implicit  Pressure  Solution. 


For  problems  involving  fluid  motion  in  spinning  cylinders,  it  is 
sometimes  convenient  to  use  stretched  meshes  that  provide  fine  zoning  near  walls 
to  resolve  viscous  shear  stresses.  In  such  cases  the  mesh  will  often  have  some 
cells  with  large  aspect  ratios  (i.e.  sx  »  $y  or  sy  >>  sx).  When  this  happens 
the  SOR  pressure  solution  technique  used  in  S0LA-V0F  may  exhibit  extremely  slow 
convergence.  To  eliminate  this  problem,  we  have  devised  a  new  solution  method 
based  on  the  Alternating  Direction  Implicit  (ADI)  scheme.  In  the  ADI  method 
the  pressure  on  a  given  row  of  cells  is  solved  implicitly  using  a  tridiagonal 
solver.  The  pressures  on  adjacent  rows  are  treated  as  known  values  and  are  not 
changed  during  this  solution  step.  Each  row  is  treated  in  this  way  in  a  cyclic 
manner.  Then  the  columns  are  swept,  one-oy-one,  using  a  tridiagonal  solver  for 
implicitly  calculating  the  pressures  in  each  column.  Switching  back  and  forth 
between  rows  and  columns  gives  the  method  its  "Alternating  Direction"  name.  If 
only  rows  or  only  columns  are  implicitly  solved,  the  method  is  unstable. 
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Our  variation  of  the  ADI  method  is  called  SADI.  In  this  scheme  we 
all ci-  some  coupling  between  rows  or  columns  during  the  one-dimensional  implicit 
solutions.  For  example,  when  solving  for  a  row  of  pressures  in  the  x-direction, 
one  half  of  the  terms  representing  pressure  influences  in  the  y-direction  are 
included.  Of  course,  pressures  in  neighboring  rows  are  still  held  fixed.  The 
"one  half"  prescription  has  empirically  been  found  to  significantly  increase 
convergence  of  the  method.  Also,  this  scheme  is  stable  when  only  rows  or  only 
columns  are  solved  for  implicitly.  The  user  selects  the  desired  pressure 
solution  method  through  the  input,  parameters  IADIX  and  I  ADI  Y.  When  I  ADI  X  =  1, 
rows  are  treated  implicitly;  otherwise  IADIX  should  be  zero.  Similarly  for 
columns,  a  value  for  I ADI Y  *  1  gives  an  implicit  solution  for  columns,  and 
I AD  I Y  =  0  otherwise. 

As  a  general  rule,  the  pressure  solution  should  be  implicit  in  the 
x-direction  if  there  are  cells  with  fix  <<  5y  and  implicit  in  the  y- 
direction  if  there  are  cells  with  fiy  <<  fix. 

3.3  Miscellaneous  Code  Changes. 

A  calculation  of  the  volume  of  "F-fluid"  in  the  computing  mesh,  FVOLT, 
has  been  added,  and  this  is  printed  out  with  each  cycle's  short  print.  F-fluid 
is  the  only  fluid  in  a  free  surface  problem,  while  in  a  two-fluid  problem  it  is 
the  fluid  with  density  RHOF.  FVOLT  is  also  printed  at  the  top  of  each  long 
print.  This  bookkeeping  feature  makes  it  easy  to  plot  the  time  history  of 
fluid  volume,  leaving  a  cylinder  with  an  open  end. 

The  initial  setup  routine  has  been  changed  so  that  fluid  can  be 
initialized  with  both  vertical  and  horizontal  interfaces.  Two  input  parameters, 
FLHT  and  FLHTV ,  control  the  fluid  distribution.  Also,  input  parameters  have 
been  added  to  assign  initial  rotation  rates  to  the  fluid  and  tc  the  walls  of  a 
cylindrical  container.  These  rates,  RPSI  and  RPSB  respectively,  are  in  revolutions 
per  second.  All  these  parameters  are  described  fully  in  Section  4. 

Some  corrections  have  been  made  to  the  PETACAl.  routine.  The  auantity 
PETA  (I,J)  is  used  as  a  pressure  relaxation  factor  for  full  fluid  cells  that 
are  interpolation  neighbors  for  surface  cells.  An  early  update  recommended  for 
this  factor  was  not  in  the  CSL  version  of  the  SOLA-VOF  code.  We  have  added  this 
correction  together  with  a  recent  modification  discovered  during  the  course  of 
this  work.  These  changes  should  solve  the  occasional  lack  of  convergence  problem 
encountered  in  the  old  code  using  the  SOR  pressure  solution  method.  Also,  a 
new  calculation  for  fr^e  surface  slopes  has  been  included,  which  works  better 
near  boundaries  and  obstacles. 

In  a  few  isolated  cases,  the  P-advection  and  NF  flag  setting  routines 
could  have  occasionally  produced  spurious  results.  Although  these  errors  are 
not  likely  to  cause  significant  effects  in  most  calculations,  we  have  corrected 
them  in  the  new  code. 

4.  USER  INSTRUCTIONS 

SOLA-VOF/CSL  is  operated  in  all  respects  like  its  predecessor  except 
for  the  addition  of  several  new  input  parameters.  These  parameters  are  described 
in  the  following  subsections.  A  complete  list  of  input  parameters  is  given  in 
the  Appendix. 


4.1 


FLHT  and  FLHTV. 


The  initial  fluid  configuration  is  defined  by  these  values  as  indicated 
in  Figure  2.  For  one  material  problem  the  fluid  is  defined  by  the  F  =  1.0 
region 


F  = 1  REGION 


Figure  2.  Initial  Fluid  Configuration 

extending  in  the  radial  or  x-direction  from  FLHTV  to  the  outer  edge  of  the 
computing  mesh  and  from  the  bottom  of  the  computing  mesh  up  to  height  FLHT.  In 
two  material  problems  the  F  =  1.0  region  has  density  RHOF,  while  the  F  =  0.0 
region  has  density  RHOFC.  Any  other  initial  fluid  configurations  must  be 
programmed  in  the  subroutine  SETUP.  If  values  for  these  parameters  are  not 
input,  they  default  to  FLHT  =0.0  and  FLHTV  =  0.0. 

4.2  AZW. 

A  parameter,  AZW,  is  used  to  indicate  when  azimuthal  velocities  are 
to  be  computed.  If  AZW  =  0.0,  the  code  will  run  in  its  original  mode  and  the 
azimuthal  velocity  equation  will  be  skipped  in  the  calculations.  To  select 
azimuthal  velocities,  set  AZW  =  1.0.  Values  of  AZW  other  than  0.0  or  1.0  have 
no  meaning.  The  default  value  for  AZW  is  0.0. 

It  should  be  noted  that  when  azimuthal  velocities  are  wanted  (i.e. , 

AZW  =  1.0),  the  geometry  must  also  be  defined  as  cylindrical,  ICYL  =  1. 

4.3  RPSB  and  RPSI. 

Initial  fluid  and/or  container  rotation  rates  can  be  specified  through 
the  parameters  RPSB  and  RPSI.  The  container  or  boundary  rotation  rate  in 
revolutions  per  second  is  RPSB.  The  fluid  is  given  an  initial  rigid  body 
rotation  in  revolutions  per  second  equal  to  RPSI.  If  there  are  obstacles 
defined  within  the  mesh  (cells  with  BETA  =  -1.0),  these  will  receive  the  rota¬ 
tion  RPSI,  which  will  then  remain  constant  in  the  obstacles  throughout  a  calcula¬ 
tion.  Boundary  rates  RPSB  will  also  remain  constant  unless  the  user  redefines 
them  within  the  code.  When  they  are  not  input,  these  parameters  default  to 
zero. 


4.4 


I AD IX  and  I ADI Y. 


To  select  the  original  SOR  pressure  solution  scheme,  the  values 
I  AD  IX  =  I AD  I Y  =  0  should  be  input  to  the  code  (these  values  are  also  the  default 
values).  Implicit  calculations  for  pressure  will  be  performed  in  the  x-direction 
if  I AD  I X  =  1  and  in  the  y-direction  if  I  ADI Y  =  1.  Both  parameters  may  be  set 
to  unity  for  the  full  SADI-type  of  alternating  direction  method.  If  either 
parameter  is  different  from  zero,  the  over-relaxation  parameter  OMG  will  have 
no  effect  as  it  is  not  used  in  SADI. 

Caution:  the  SADI  scheme  will  not  work  with  periodic  boundary  condi¬ 
tions.  Periodic  boundaries  require  a  more  complex  solution  strategy,  which 
has  not  been  implemented  in  SOLA-VOF/CSL.  SADI,  however,  will  work  with  any 
combination  of  free  boundaries  and  interior  obstacles. 

4.5  INTDC. 

This  Parameter  controls  the  difference  approximation  used  in  advecting 
the  F-function  that  defines  free  surfaces  and  fluid  interfaces.  In  free  surface 
applications,  the  advection  scheme  should  use  an  acceptor-cell  type  of  approxima¬ 
tion  in  the  interior  of  F  =  1.0  regions.  This  is  insured  by  setting  INTDC  =  0. 

For  two-fluid  applications,  it  is  usually  best  to  use  a  donor-cell 
approximation  in  the  interior  F  regions.  This  is  obtained  by  setting  INTDC  =  1. 

The  correct  settings  for  INTDC  are  automatically  established  within  the  code  as  a 
default.  Thus,  INTDC  should  only  be  input  when  the  user  wishes  to  override  the 
default  values. 

5.  SAMPLE  PROBLEMS 

5.1  Free  Surface  Example. 

As  an  illustration  of  the  new  code's  capabilities  for  computing  the 
flow  of  fluid  leaving  a  spinning  cylindrical  container,  we  have  run  the  following 
problem.  A  cylinder  of  radius  R  =  1.0  and  length  L  =  3.3  was  initialized  with 
fluid  occupying  80%  of  its  volume  (20%  void).  Nondimensional  units  have  been 
used  for  convenience.  The  fluid  and  cylinder  are  initially  in  solid  body  rotation 
at  92  revolutions  per  time  unit.  No  axial  gravity  was  assumed  so  that  the 
fluid  has  a  free  surface  parallel  to  the  axis  of  the  cylinder.  At  time  zero 
the  top  of  the  cylinder  has  been  instantaneously  raised  by  0.7  distance  units. 

A  constant  pressure  boundary  condition  is  assumed  at  the  outer  radius  in  the 
gap  formed  when  the  top  was  raised.  Figure  3  shows  the  mesh  setup  used  for 
this  problem.  Figure  4  contains  several  velocity  plots  obtained  in  the  course 
of  the  calculation.  A  plot  of  the  volume  of  fluid  contained  in  the  cylinder 
versus  time  is  given  in  Figure  5.  This  problem  took  a  little  less  than  10 
minutes  on  a  CDC-750  computer. 

For  new  users  the  following  discussion  will  clarify  how  the  code  was 
set  up  to  perform  this  calculation.  The  finite-difference  mesh  must  be  constructed 
for  the  cylinder  with  raised  top  and  to  have  only  a  portion  of  the  r  =  R  boundary 
open  for  outflow.  The  easiest  way  is  to  set  the  right-hand  boundary  condition 
as  a  specified  pressure  boundary  and  to  define  obstacle  cells  in  the  right-most 
column  of  cells  over  that  portion  of  the  boundary  that  represents  the  cylinder 
wall.  The  remaining  cells  in  that  column  (above  the  obstacle)  will  then  have 
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Figure  3.  Mesh  Setup  for  Example  Problems 
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Figure  4.  Velocity  Fields  for  Example  Problem  A 


Time  (T)  and  maximum  velocity  (VM)  in  each  frame  is: 
T  =  .004,  VM  =  293:  T  *  .016,  VM  =  262:  T  =  .089,  VM 


I.M  2. BO  3.00  4  00  S.M  B.M  7.00  8.00  0  00  10.00 
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igure  5-  Volume  of  Fluid  in  Cylinder  Versus  Time 
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the  specified  outflow  boundary  pressure,  p  «  0.  To  account  for  this  obstacle 
we  must  define  the  outer  mesh  radius  to  extend  one  mesh  cell  further  than  the 
radius  of  the  cylinder.  Using  the  automatic  mesh  generator  in  the  code,  we  can 
easily  do  this  by  inputting 

NKX  =  1 

XL ( 1 )  =  0.0  XL( 2 )  =  1.05 

XC(1)  =  1.0 

NXL(l)  =  7 

NXR(l)  =  1 

OXMN(l)  =  0.05 

This  input  assigns  XC(1)  as  the  cylinder  radius  and  Xl(2)  as  the  radius  plus 
one  cell,  NXR(l)  =  1,  with  widtn  DXMN(l)  =  0.05.  In  the  sample  problem  shown 
in  Figure  4,  we  wanted  to  have  several  cells  of  size  0.05  adjacent  to  the 
cylinder  wall  so  we  decreased  XC(1)  and  increased  the  number  of  cells  between 
XC(1)  and  XL(2), 

XC(1)  =  0.85 
NXR(l)  =  4 

For  the  axial  direction,  we  wanted  small  cells  at  the  bottom  of  the  cylinder 
with  uniform  stretching  toward  the  top.  Thus,  we  set 

NKY  =1 

YL(1)  =  0.0  YL(2)  =  4.0 

YC  ( 1 )  =  0.5 

NYL(l)  =  5 

NYR(l)  =  25 

DYMN(l)  =  0.1 

which  puts  5  cells  with  the  minimum  height  0.1  at  the  bottom,  and  30  cells 
altogether  to  extend  the  mesh  to  a  maximum  height  of  4.0.  In  this  problem  the 
height  corresponds  to  the  distance  between  the  cylinder's  bottom  end  its  top  in 
the  raised  position.  We  selected  3.3  as  the  initial  height  of  the  cylinder 
because  this  corresponded  to  the  top  of  cell  j  =  27  in  the  mesh.  Thus,  we 
defined  the  cylinder  wall  by  setting  the  last  column  of  cells,  i  =  12,  from  j  =  2 
to  j  =  27  as  obstacle  cells,  i.e. ,  we  set  BETA  (i,j)  =  -1.0  in  these  cells. 

This  is  done  in  the  SETUP  routine. 


More  generally,  to  define  the  vertical  mesh  so  that  it  will 
have  a  mesh  line  at  the  top  of  the  cylinder  wall,  y  =  y^y,  and  the  raised  top  at 
y  =  yy,  the  y-mesh  input  should  have  the  form: 


NKY  =  2 

YL(1)  =  0.0,  YL(2)  =  ycw,  YL(3) 
YC(1)  =■  yB,  YC(2)  5  yT  -  yT 
NYL(l)  *  1,  NYL(2)  =  ? 

NYR(l)  =  ?,  NYR(2)  =  1 
DYMN(l)  =  yB,  DYMN(2)  =  yT 


=  n 


where  yB  is  the  cell  size  wanted  at  the  bottom  wall  of  the  cylinder  and  yy  is 
the  cell  size  wanted  at  the  top.  The  cell  numbers  NYL(2)  and  NYR(l)  depend  on 
the  choices  for  yg,  yy,  and  the  overall  mesh  resolution  that  the  user  wants. 
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To  initialize  the  fluid  in  the  mesh  for  an  80%  fill,  we  first  input 
FLHT  =  3.3  so  that  the  fluid  extends  up  to  the  top  of  the  closed  cylinder. 

Then  FLHTV,  which  controls  the  location  of  the  left  side  of  the  fluid  region, 
must  be  calculated  to  give  the  80%  fill  (or  20%  void),  i.e. ,  FLHTV  =  0.20R2, 
where  R  is  the  cylinder  radius.  In  the  present  example,  FLHTV  =  0.4472. 

Because  we  assume  an  initial  solid  body  rotation  the  boundary  and 
fluid  rotation  rates  are  equal,  RPSI  =  RPSB  =  92. 

For  problems  with  rotation  we  must  also  input  AZW  =  1.0  and  ICYL  -■  1. 
The  boundary  conditions  must  have  WL  =  1  because  the  left  boundary  is  the  axis 
of  rotation,  WT  =  WB  =  2.0  to  have  no-slip  conditions  at  the  cylinder  ends,  and 
WR  =  5,  which  gives  the  specified  pressure  condition  along  that  portion  of  the 
rignt  boundary  that  is  open  for  flow-  The  obstacle  along  the  lower  portion  of 
the  right  boundary  defines  the  cylinder  wall  and  keeps  the  initial  solid  body 
rotation  rate  assigned  to  it. 

All  the  remaining  input  data  required  for  this  problem  consist  of  the 
controls  for  frequency  of  output  data,  time  limits,  and  initial  time  step.  For 
the  initial  time  step  we  used  0.0001  as  a  guess;  however,  the  code  cuts  this 

down  somewhat  as  fluid  begins  to  move.  We  also  selected  a  viscosity  of  NU  = 

0.05.  This  choice  was  arbitrary  and  was  used  only  to  have  non-zero  viscous 
stresses  in  order  to  check  the  code  with  this  option. 

5. 2  Two-Fluid  Example. 

As  an  illustration  of  a  rotating  system  with  two  fluids,  we  have  used 
nearly  the  same  setup  as  for  Example  A.  Because  free  surfaces  are  not  allowed 
when  using  the  two-fluid  option,  the  top  of  the  cylinder  has  been  moved  down  to 
a  closed  position.  This  is  done  in  the  input  by  defining  YL(2)  =  3.3  and 
reducing  the  cell  number  NYR(l)  from  25  to  21.  For  consistency  we  have  retained 

the  obstacle  cells  in  column  i  =  12  and  WR  =  5.  It  would  be  better,  in  general, 

to  eliminate  the  obstacle  by  reducing  the  mesh  in  the  x-direction  to  extend 
only  out  to  the  cylinder  wall,  and  using  a  WR  =  2  boundary  condition  at  the 
right  boundary. 

The  initial  fluid  configuration  consists  of  a  fluid  with  density  1.0 
(RH0F  =  1.0)  filling  the  lower  half  of  the  cylinder  (FLHT  =  1.65,  FLHTV  =  0.0) 
and  a  fluid  of  density  0.9  (RH0FC  =  0.9)  filling  the  top  half.  The  system  is 
initially  in  solid  body  rotation,  RPSI  =  RPSB  =  92.  Thus,  we  envision  the 
fluids  having  been  spun  up  to  a  steady  state  rotation,  which  requires  some  sort 
of  diaphragm  to  keep  them  from  mixing.  At  the  start  of  our  calculation,  this 
diaphragm  has  been  removed  and  we  wish  to  see  how  the  two  fluids  will  begin  to 
mix. 


Figures  6  and  7  show  the  early  fluid  interface  deformation  and  velocity 
fields  that  were  calculated.  Initially  the  heavier  fluid  on  the  bottom  wants 
to  run  up  the  outer  wall  into  the  lighter  fluid.  This  trend  only  persists  for 
a  short  time,  however,  because  the  large  eddy  flow  associated  with  this  displace¬ 
ment  tries  to  carry  fluid  near  the  outer  wall  inward,  and  fluid  in  the  inner 
region  outward.  Angular  momentum  carried  with  this  exchange  produces  an  unstable 
situation  and  results  in  complicated,  and  highly  transient,  secondary  flows.  At 
the  time  of  the  last  plot,  in  fact,  the  inner  region  has  undergone  a  complete 
flow  reversal.  The  initial  counterclockwise  flow  has  reversed  to  a  large 
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Figure  7.  Velocity  Field  for  Example  Problem  B 

Time  (T)  and  maximum  velocity  (VM)  in  each  frame  is:  T  =  0.0, 
VM  =  0.0;  T  =  .005,  VM  =  31 ;  T  =  .01,  VM  =  26;  T  =  .015,  VM  =  37. 


21 


clockwise  eddy  extending  out  to  approximately  0. 9R  with  a  long,  narrow  counter¬ 
clockwise  flow  along  the  top  portion  of  the  outer  wall.  With  this  flow  reversal, 
the  fluid  interface  near  the  axis  is  now  moving  upward. 

Approximately  7  minutes  of  CDC-750  computer  time  was  required  to 
obtain  these  results.  The  problem  Is  relatively  slow  running,  probably  because 
of  its  highly  transient  character.  Lack  of  funds  in  our  computer  budget  prevented 
us  from  running  this  problem  further.  However,  it  has  certainly  gone  far 
enough  to  illustrate  this  new  capability  of  S0LA-V0F/CSL. 
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APPENDIX 


SOLA-VOF/CSL  INPUT  DATA 


The  following  list  describes  the  input  data  for  SOLA-VOF/CSL.  Default 
values  are  indicated  for  all  quantities.  Default  values***  indicate  quantities 
that  must  be  specified  for  each  calculation. 


NAMELIST/XPUT 


Quantity 

Default 

DESCRIPTION 

ALPHA 

1.0 

Controls  donor/centered  fluxing 

AUTOT 

1.0 

Selects  automatic  time  step  control 

AZW 

0.0 

Set  to  1.0  for  an  azimuthal  velocity 

CANGLE 

90 

Contact  angle  in  degrees 

CSQ 

-1.0 

Material  soundspeed  (-1.0=incompressible  flow) 

DELT 

Time  step  size 

EPSI 

.001 

Pressure  convergence  criterion 

FLHT 

0.0 

Fluid  surface  in  y-direction 

FLHTV 

0.0 

Fluid  surface  in  x-direction 

GX 

0.0 

Body  acceleration  in  x-direction 

GY 

0.0 

Body  acceleration  in  y-direction 

I  AD  IX 

0 

Set  to  1  for  SADI  in  x-direction 

IADIY 

0 

Set  to  1  for  SADI  in  y-direction 

ICYL 

0 

Set  to  1  for  cylindrical  coordinates 

IMOVY 

0 

Set  to  1  for  movie  output 

Set  to  1  for  interior  F  donor  cell 
advection 


INTDC 


0 


Default 


Description 


i 

|  ISURFIO 

0 

Set  to  1  for  surface  tension 

X 

ISYMPLT 

0 

Set  to  1  for  symmetric  plot 

NMAT 

1 

Set  to  2  for  two-fluid  applications 

1 

NPX 

0 

Number  of  particles  in  x-direction 

NPY 

0 

Number  of  particles  in  y-dlrectlon 

NU 

0 

Kinematic  viscosity 

1 

OMG 

1.7 

Over-relaxation  factor 

.V 

PLTDT 

Time  between  plots 

PRTDT 

Time  between  prints 

3 

RHOF 

1.0 

Density  of  F=1  fluid 

*  ■ 

RHOFC 

1.0 

Density  of  F=0  fluid 

RPSB 

0.0 

Revolutions  per  sec  of  boundary 

i 

RPSI 

0.0 

Revolutions  per  sec  of  fluid  (initially) 

X 

•> 

SIGMA 

0.0 

Surface  tension  coefficient 

•  w* 

i 

i 

TWFIN 

Problem  finish  time 

UI 

0.0 

Initial  u-velocity 

„  J 

VI 

0.0 

Initial  v-velocity 

V 

V 

'i 

VELMX 

1.0 

Scale  for  velocity  plots 

WB 

1 

Bottom  boundary  condition 

*v 

W'. 

1 

Left  boundary  condition 

WR 

1 

Right  boundary  condition 

« 

WT 

1 

Top  boundary  condition 

O 

XPL 

0.0 

Left  side  of  particle  region 
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Quantity  Default 


Description 


XPR 

0.0 

Right  side  of  particle  region 

YPB 

0.0 

Bottom  side  of  particle  region 

YPT 

0.0 

Top  side  of  particle  region 

NAMELIST/MSHSET 

DXMN(N) 

★  ★★ 

Minimum  cell  size  in  x-direction 

DYMN(N) 

Minimum  cell  size  in  y-direction 

NKX 

★  ★★ 

Number  of  x  submeshes 

NKY 

*★* 

Number  of  y  submeshes 

NXL(N) 

★  ★★ 

Cells  between  XL(N)  and  XC(?«) 

NXR(N) 

*★* 

Cells  between  XC{N)  and  XL(N+1) 

NYL(N) 

★  ★★ 

Cells  between  YL(N)  and  YC(N) 

NYR(N) 

Cells  between  YC(N)  and  YL(N+1) 

XC(N) 

★  ★★ 

x  coordinate  of  smallest  cell  location 

XL  ( N ) 

★  ** 

Left  side  of  submesh  N 

YC(N) 

★  ★★ 

y  coordinate  of  smallest  cell  location 

YL(N) 

**★ 

Bottom  of  submesh  N 
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